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1. Introduction 


Quasi-elliptic schemes arise, for example, when central differencing is used to 
approximate odd-order derivatives m elliptic systems of partial differential equa- 
tions, such as the Cauchy-Riemann, Stokes and Navier-Stokes systems Usual 
finite element approximations to such systems also lead to quasi-elliptic schemes. 
Such schemes are in some sense unstable certain highly-oscillating components are 
amplified in the discretized solution much more than in the differential solution 

Instead of the quasi-elliptic schemes, other discretizations of the same system 
can usually be constructed which are h-elliptic, hence fully stable, and which are 
also more accurate than the quasi-elliptic schemes Sometimes, however, these 
fully elliptic schemes are inconvenient to use In case of elliptic systems with odd- 
order derivatives, for example, full ellipticity is obtained by grid staggering. 1 e., 
by approximating different functions on different grids (cf [3] and [8]) This is 
inconvenient, especially near curved boundaries. Also the instability of quasi- 
elliptic approximations seldom really hurts, since the unstable components have 
very small amplitudes, which are still small even in the discrete solution The 
inaccuracy is modest The error in the quasi-elliptic solution is typically twice to 
four times larger than the error m an elliptic solution using the same grid size 
Thus, quagi-elliptic schemes are often preferred and are widely used 

The instability of quasi-elliptic schemes does seem to hurt when multigrid 
solvers are applied - The asymptotic convergence turns out to be slow, and a 
simple mode analysis traces this slowness to the unstable modes. One approach, 
perhaps the best, to deal with this difficulty is simply to ignore it - the algebraic 
slowness does not matter because it occurs m modes whose amplitudes in the 
algebraic solution are erroneous anyway, bearing no relation to their amplitudes 
in the true differential solution One should only take care not to initially admit 
large unstable amplitudes, and to average them out m case they must latter enter 
We show, by mode analyses and numerical experiments, that the usual FMG 
algorithm is very effective in solving quasi-elliptic problems to truncation level 
(i.e., to the point where algebraic errors are dominated by discretization errors) 
Sometimes the FMG solution may even be better than the exact solution of the 
discrete equations, because the unstable components of the latter are slow to enter 

Although this is the easiest approach for obtaining fast differential conver- 
gence (convergence to the differential solution through a sequence of grids), an- 
other algorithm is presented below which does provide fast algebraic convergence 
for quasi-elliptic schemes This algorithm, based on multiple coarse-grid correc- 
tions, is interesting in its own right, since it is the simplest example of a new 
kind of algorit hms for solving problems with highly-oscillating solutions, including 
highly indefinite problems (see [1, §3.2], [8] and a subsequent article) Smoothing 
rate analysis, for one quasi-elliptic example, suitably modified to account for the 
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multiple coarse-grid corrections shows that the new algorithm should algebraically 
be as efficient as usual multigrid cycles are for fully elliptic schemes Numerical 
experiments exactly yield these expected convergence rates (see §8 6, tests of such 
algorithms were also reported in [6] ) 

The significance of the present studies goes beyond elliptic PDE systems 
many non-elliptic systems, such as all subsonic steady-state flow problems, have 
determinants with at least one elliptic factor Most discretizations of such sys- 
tems provide quasi-elliptic approximation to that factor, leading to troubles and 
requiring cures similar to those reported here 

Moreover, the techniques described m this article illustrate the following gen- 
eral multigrid approaches to general non-elliptic problems (l) Differential not 
algebraic convergence is sought and usually easily obtained Modified methods 
for apnori analyzing and aposteriori measuring such a convergence have been de- 
veloped (n) With considerably more effort, fast algebraic convergence can also be 
obtained (in) The analysis of difference schemes, and the derivation of efficient 
smoothers, for any PDE system is based on the factors of the h-prmcipal part of 
the operator determinant 

We thank Ruth Golubev for some of the calculations reported m §8 


2. Definitions and Examples 


In the following L h will represent a system of q real difference operators on 
q grid functions where h, the meshsize of the grid, is for simplicity assumed 
to be uniform and the same m all directions That is, L h is a q x q matrix of 
real polynomials in T \ , , , T^ 1 , where T t are the grid translation 

operators defined by 

T\' Tj d u(x) = u{x 4- vh), 

with x = (zj , , Xd), v = (t'l , , I'd) and d being the dimension of the Euclidean 

space housing the grid (In case of staggered grids there may appear non-mtegral 
powers of T 3 and L h will most usually be a matrix of polynomials m T 1 ’ 2 and 

j = i ,<o 

Three common examples of difference operator are 
(i) The five-point (compact) Laplacian 


— -^{T 0 i -I- r lt0 + To,_i -t- T-1,0 


4To,o) 



1 ' 
-4 1 , 

1 


(2 1 ) 


where T a ^ = T°T® and the array on the left is the usual pictorial 


description 
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of the weights of the operator This is the simplest approximation to the two- 
dimensional Laplace operator A'= d 2 /dx\ + d 2 fdx\ 

(ii) The central non-staggered approximation to the Cauchy-Riemann opera- 
tor 



where d‘ = i(T, - T" 1 ) 

(iii) The central non-staggered approximation to the Stokes operator in two 
dimensions 



For simplicity we will deal m this article only with constant-coefficient oper- 
ators L h In this case the symbol L h (0) of L h is defined by 

L h Ae x - * /h = L h {6)Ae'? (|0| < 7 r) 

for any g-vector A, where 0 = (0 1} 0 x = 0iXj + • • • + and |0| = 

max(|0i|,. .,|^d|) Thus, L h (0) is a g x g matrix of polynomials in e ±x6] , j = 
1, . , d, obtained from L h by replacing each T } with . 

Also for simplicity we will deal here only with homogeneous operators L h , 
1 e , operators for which all terms in det L h (the determinant of L h ) have the 
same power in h (This means that L h approximates a homogeneous differential 
operator L, 1 e., det L is a homogeneous polynomial in d/di\,.. ,d/dxd All 
examples above are homogeneous) For homogeneous difference operators, the 
general notion of ellipticity measure on a given scale (cf [2, Sec. 3 lj or [3, Sec. 
2.1]) is not needed, and we can use the following simpler definition. 

Definition The homogeneous difference operator L h is elliptic of order 2m 
iff 

d 

I det i h {0) |> C/r 2m for a11 l£ l< (2 4) 

3 = 1 

where C is positive and independent of 0. 

Ellipticity of differential operators is defined m the same way (The parameter 
h is arbitrary then, and the range of 0 is unrestricted. It is thus more natural in the 
continuous case to replace 6/h by another phase variable, u = 0/h say.) It is easy 
to see that both A and A h are second-order elliptic. Generally, simplest central 
approximations to second-order scalar (g = 1) elliptic operators axe themselves 
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elliptic. But not all central approximations are For example, the ‘‘skew Laplacian” 

1 0 r 

0-4 0 
.1 0 1. 

or the “long Laplacian” 

= ^2(72.0 + To, 2 + T-2,0 + To, -2 - 4 To, 0) (2 6 ) 

both approximating A, have the symbols 

A* (0) = [(cos 6 1 — cos 62 ) 2 + sin 2 9\ + sm 2 62} 

A 2h {0) = -^-(sin 2 9\ -f sm 2 B2) 
h z 

which clearly fail to satisfy (2 4) Indeed. A x (x. n) =0 and A 2h (7t. 0) = A 2/l (0,7r) 
= A 2h (n, 7r) = 0 Whereas these examples seem somewhat artificial (although the 
skew Laplacian does naturally arise in various situations, eg., in semi-implicit 
Lagrange codes [4, § IV] and for some kinds of finite elements [7]) non-elliptic 
operators are very common m approximations to elliptic systems ( q > 1) The 
discrete Cauchy-Riemann (2.2) and Stokes (2 3) operators well represent this sit- 
uation- They are the simplest (non-staggered) central approximations to elliptic 
operators, but det Lq R = A 2h and det Lg — A h A 2h , hence they do not satisfy 
(2 4), their symbol vanishing wherever A 2h does Note that taking the determi- 
nant commutes with passing to the symbol, hence ellipticity of L h is equivalent to 
ellipticity of det L h , which in turn is equivalent to ellipticity of all factors of det 
L h . 

Finite element discretizations of the same elliptic systems, with uniform non- 
staggered partitions, give rise to similarly non-elliptic difference operators This is 
not usually recognized because finite element discretizations are seldom Fourier- 
analyzed as uniform-grid operators 

In all the above examples, even when L h fails to satisfy (2 4), it still satisfies 
the weaker condition 

d 

| det L h (9) |> Ch~ 2m Y^ sm 2m (^), for all \6\<n, (2 7) 

i = i 

where C is positive and independent of 0 The term quasi- elliptic was introduced 
in [8] to describe such operators 

Perhaps all reasonable approximations to homogeneous elliptic equations sat- 
isfy (2.7), but for the purpose of including some additional, not-so-reasonable ap- 
proximations, we can extend the class of operators, and admit any homogeneous 



1,1 + Ti-! 4- T- 1,1 + T_ i,— i - 4T 0 ,o) = ^ 


2h? 
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operators L h for which det L h (6 ) vanishes only at a finite number of points. This 
class includes for example A 2hx = (72,2 +72,-2 + 71-2 , 2 + T- 2,-2 - 47b,o)/(8h 2 ), 
which satisfies neither (2 4) nor (2.7), but for which the methods described below 
are still applicable 

More generally when inhomogeneous operators are also admitted, our meth- 
ods will extend to any operator L h with 0(1) “measure of quasi-ellipticity” , defined 
by 

E h ’°{L h )= min | det L h {6) | /| det L h {6') |, (2.8) 

a>|0|>|0'l 

for some reasonable a > 0 E h ’ n is the usual measure of elhpticity E h described 
in [3] The methods here w'ill in principle work for any positive a, although they 
will gradually deteriorate with the decrease of a for which E h,a (L h ) is still 0(1). 

For clarity, we discuss below only homogeneous operators, and the strict quasi- 
ellipticity (2 7) is assumed. 


3. Instability and Inaccuracy 

Quasi-elliptic operators do meet some general stability requirements even if 
they do not satisfy (2 4) For example, the skew Laplacian (2 5) is a positive type 
operator, hence satisfying the maximum principle. The associated matrix has a 
dominant diagonal Nevertheless, m a certain sense such operators are not quite 
stable. Namely, since det L h (0) = 0 for some 0 + 0, in an infinite space, or under 
periodic boundary conditions, there exists a highly-oscillating function v h (x) = 
Aexp(i9 x/h) which satisfies the homogeneous equation L h v h (x) = 0. Hence 
the solution, unlike the corresponding differential solution, is not unique (upto 
an additive constant); it contains an undetermined highly-oscillating component. 
Similarly, in any bounded domain with any boundary conditions, functions w h (x) 
close to v h (x) (e.g . w h = <P\v h + tp?, being smooth) exist which satisfy the 
boundary conditions and for which L h w h is everywhere very small. Such w h 
therefore forms an unstable mode. A small change in the equation can introduce 
a large change proportional to w h . This is a kind of numerical instability, since a 
corresponding large change in the differential solution cannot occur 

This numerical instability need not hurt much - If the differential system is 
LU = F and the discrete system is L h U h = F h , all one has to do is to define 
F h — I h F, say, through an averaging operator I h which liquidates the unstable 
modes, i.e I h {0) = 1 and the ratio I h {Q)/L h (6) is uniformly bounded for all 
|0| > e > 0 For example, one can take I h = S h I ,h , where I' h is any F averaging 
suitable for the fully elliptic case and S h is like the solution averaging S h described 
below. Even this is unnecessary in the usual, smooth case (m the same way that 
the above rule for I h is frequently neglected for fully elliptic L h ), because the 
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unstable modes, even when unduly magnified by the discretization, are usually 
still small 

Generally, the main disadvantage of quasi-elliptic operators is a certain loss of 
accuracy compared to corresponding truely elliptic operators, which is simply due 
to the larger differencing steps taken in certain terms of the quasi-elliptic scheme 
In some cases this is particularly obvious, since the grid is locally decoupled into 
several subgrids which are not connected to each other by the quasi-elliptic oper- 
ator For example, the skew-Laplacian (2.5) introduces no coupling between red 
and black points (in the usual sense of checkerboard coloring, one color being as- 
sociated with gridpomts where (xi + )/h is odd, the other with even) On each 

subgrid the discretization looks like the compact Laplacian (2 1) on a rotated grid 
with meshsize h\ = y/2 h Similarly, in case of (2 2). the grid is decoupled into 
4 staggered subgrids with meshsizes 2h, on each of which the operator has good 
ellipticity (being m fact equivalent to the staggered-grid approximation described 
m [3, §17.2] or [5. §5 2]) Thus, since the approximation is 0(h 2 ). the error in 
case of (2 5) is on the average twice larger, and in case of (2 2) four times larger 
than the errors m corresponding fully elliptic approximations (assuming other dis- 
cretization errors, related for example to the representation of right-hand sides or 
boundary conditions, behave similarly) In these cases, in other words, each of 
the subgrids can produce the resulting accuracy by itself, other subgrids only add 
work 

When derivatives are calculated from the solution, however, the approximat- 
ing difference quotients may show much greater loss of accuracy, because they 
involve differences between values belonging to different subgrids The error m 
f-order derivatives will generally be 0{h~ l ) times the errors in the function itself 
This excessive error can be avoided by taking differences only from one subgrid at 
a time, or, more generally, by using only difference operator D h such that D h {6) 
vanishes wherever L h {6) does In case of (2 3), for example, derivatives of the 
third unknown function (the pressure) should be approximated by long differences 
such as d*. d c ,di. etc. 

The instability described above can also be removed, and the inaccuracy in 
derivatives proportionally reduced, by averaging the solution, that is, by replacing 
the computed solution u h by S h u h , where S h is an averaging operator which 
removes all the unstable components In other words, S h (0) = 1 and outside 
a neighborhood of 6 = 0 the ratio S h (6)/L h (0) should be uniformly bounded 
(wherever defined) For the quasi-elliptic L h satisfying (2.7) there always exists 
such an averaging operator of the form 

s h = nj =1 + |r/ 1/2 ) mj , (3.i) 

with integral m } < m. On the other hand, the averaging may further reduce the 
accuracy of the solution. With the averaging (3 1) the lost accuracy is 0(h 2 ) One 
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can make that loss 0(h 2a ) by taking for example 

S h = nj =1 (l - (i - i T, - (3 2) 

Another slight difficulty typical to quasi-elliptic approximations is the need to 
define extra boundary relations This can satisfactorily be done by extrapolation 
(cf , e.g., §8 2) 

In summary, although quasi-elliptic discretizations are in principle inferior 
to fully elliptic ones (obtainable for systems by grid staggering), they can be 
used Since many programmers consider grid staggering a serious complication, 
especially near general boundaries, quasi-elliptic schemes and their fast solution 
become important 


4. Multigrid Troubles and Their Implications 


Usual multigrid solvers yield poor asymptotic convergence rates when applied 
to quasi-elliptic schemes (see [4] and §8 4 below). The reason is simple. Slow to 
converge are the unstable modes, such as v h or w h above. They cannot signifi- 
cantly converge by coarse-grid corrections, since they are high-frequency modes, 
essentially invisible on coarser levels Neither can they significantly converge by 
any type of relaxation, since an error like w h shows a very small residual function 
L h w h (compared with residuals shown by other modes with comparable ampli- 
tude) and the corrections introduced by any relaxation scheme are proportional 
to the size of the residuals (cf. [3, Sec 1 1]) In particular, the amplification 
factor of the error mode exp(i# x/h) per relaxation sweep must be 1 when 
L h (6) = 0, and since the latter equality holds for some |£i = 7r, the smoothing 
factor /2 = max„/ 2 <|e|<- 7 r |m(£)I cannot be smaller than 1. 

The poor asymptotic rates are not a real trouble , though The modes slow to 
converge are exactly those unstable modes for which algebraic convergence is not 
really desired, their amplitudes in the algebraic solution being unrelated to their 
amplitudes in the differential solution The only concern is that these amplitudes 
will remain suitably small 

This situation is typical to all problems which are not fully elliptic, includ- 
ing most problems in fluid dynamics. Slow asymptotic convergence of suitable 
multigrid cycles occur exactly in those components where not much convergence 
is needed anyway Whenever this situation arises, it is m a sense an absurd to 
try and fix the algorithm (although we show in Sec. 7 below how to do it), since 
one would then often end up investing most of his human and computer resources 
to obtain improvements which are meaningless in terms of solving the original 
differential equations. 
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Thus, the real objective of.multigrid solvers should not be a fast algebraic 
convergence (convergence of the computed solution u h to the exact discrete solution 
U h ), but fast differential convergence (convergence of u h to the true differential 
solution U), using any sequence of meshsizes h and measured directly m terms 
of the decrease in || u h — U || as function of the overall computational work (cf 
[3, §13]). This modified objective allows for simpler algorithms, but also calls for 
some modifications in our approach for analyzing algorithms, for aprion predicting 
and aposterion measuring their performance The next two sections will illustrate 
these modifications for the case of quasi-elhptic schemes 


5. Modified Mode Analysis 


It was shown m Sec 4 that in case of quasi-elhptic systems fj > 1, but that 
this bad smoothing factor is not relevant to our real objective To analyze a given 
relaxation scheme, assume first that it is as efficient as needed for the differential 
convergence of the highest frequency modes (which should latter be checked by the 
2-level FMG mode analysis mentioned below) The question then is what efficiency 
one should expect from the multigrid cycle (employing the given scheme on all 
levels) in reducing all other modes As in the conventional smoothing analysis, our 
simplifying assumption here will be that relaxation on each level should efficiently 
treat all modes in only one segment of modes, and that the union of these segments 
should cover all relevant modes Instead of assigning to grid h the conventional 
segment n/2 < |0| < 7r, however, we can assign to it any segment of the form 
a/2 <|| 8 j|< a, with any norm |j 6 ||. That would automatically assign to grid 
h/2 the segment a/4 <|| 0 ||< a/2, and so on It means that we allow some of the 
highest frequency components on any intermediate level not to converge efficiently 
by relaxation on that level, as long as those components efficiently converge by 
the next-finer-level relaxation This only leaves the highest frequency modes on 
the finest grid unaccounted for, which is exactly the segment where we do not 
seek simple algebraic convergence Thus, the modified definition of the smoothing 
factor relevant for our purpose here is 

u = min max |u(0)|, (5 1) 

a/2<||e|l<a 


where fi( 0 ) is the amplification factor of exp(i0 x/h) per relaxation sweep, and 
the minimum can be taken over all a > 0 and over all possible choices of the norm 
|| • || (For a generalization of this definition to cases of semi coarsening, cf [3, 


§ 12 ]). 


In case of the skew Laplacian (2 5), for example, the lexicographically order 
Gauss-Seidel relaxation yields the amplification factor 


n( 0 ) = e%6i cos 62 / (2 — e ' e> cos 0 2 ), 


(5.2) 
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so that n(n, 7r) = 1 and the conventional smoothing factor is 1. But choosing 
a — 7r and || 0 ||= max(|0i + 02 1, |0i — 02 i) easily shows that the modified factor 
(5 1) yields p < 5. The same can be shown for the long Laplacian (2.6), by taking 
a = 7r/2 and || 0 || = |£! 

In case of systems ( q > 1). the q x q amplification matrix n of the mode 
Aexp(i0-x/h ) depends both on-0 and on the q - vector A The modified smoothing 
factor p is then defined by 

p — mm max |j pA || / || A || 

a/ 2< ||0|j <a 

where a is allowed to depend on both 0/10! and A/|A| With these definitions and 
suitable distributed Gauss Seidel (DGS) relaxation schemes (see e g [3, §18.6]) 
this again yields p < 5, for both the Cauchy-Riemann (2.2) and the Stokes (2 3) 
operators In all these cases, still better factors are obtained by four-color ordering, 
for which definitions (5 1) and (5 2) should further be extended (cf (3 2) m [3]) 

As for two-level analyses (cf [3. §4 lj or [5, §4.6]). they always couple lowest 
with highest frequency modes In non-elhptic cases some highest-frequency modes 
are not expected to converge fast What the analysis should then tell us is how 
efficient is the entire multigrid algorithm m reducing the algebraic errors below 
the truncation errors This can be done by a two-level FMG mode analysts, which 
Fourier analyzes the A T -FMG algorithm described below (usually for N = 1) by 
assuming exact solution of the coarse grid equations (both for obtaining the first 
approximation and in each of the N cycles) and by comparing for each mode the 
final algebraic error with the truncation error (see [3, §7.4]) 


6. FMG Solution to Truncation Level 

Since the multigrid cycling is inefficient m reducing unstable mode errors, 
the multigrid solver should take care not to start with an initial solution which 
contains large amplitudes of such errors The overall initial error in unstable modes 
should better be smaller than the overall truncation error. This is easily obtained 
by taking a first approximation from a coarser grid, employing interpolation of 
suitable order The usual “Full multigrid” (FMG, also c allied “nested iteration”) 
algorithm can therefore be used, with slight modifications The usual algorithm 
and its modifications are briefly described in the following For a flowchart, and 
a detailed discussion of FMG algorithms and the order of the first interpolation, 
see Secs 1.6 and 7 in [3] For simplicity we describe here the Correction Scheme 
(CS) version of the algorithm, so the problems are assumed linear; it should be 
converted to Full Approximation Scheme (FAS) to treat nonlinear problems [3, 
§ 8 ] 
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6.1 Multigrid cycle A sequence of grids is given with meshsizes h k (k = 1,2,3, 
), where h k +i = h k /2 On the h k grid the discrete equations have the form 

L k U k = F k (6 1) 

where L k approximates Z, fc+1 Given u k , an approximate solution to (6.1), the 
multigrid cycle MG for producing an improved approximation u k 

u k «- MG{k,u^F k ) (6.2) 


is recursively defined as follows 

If k = 1 solve (6 1) by any direct or iterative method, yielding the final result 
u k Otherwise do (A) through (D) 

(A) Perform u\ relaxation sweeps on (6 1), resulting in a new approximation 
u k 

(B) Starting wuth Ug~ ] =0 make 7 successive cycles 

u k - } - MG(k - l.u k :}.I k ~ 1 (F k - L k u k )). (j = l.. .7) 

where /£ -1 is a transfer (‘‘reduction”) of residuals from grid h k to grid h k ~i We 
have used the “full weighting” 

I k ~ 1 = \ T o.o - g(?b,i + T h0 + T 0 -i + T- h o) 

+ ie(Ti,i + Ti-i T_i,i + 

(C) Calculate u h = u k where /£_ j is a suitable interpolation 

(“prolongation”) from grid h k - 1 to grid h k For problems considered here, bilinear 
interpolation is used 

(D) Perform relaxation sweeps on (6 1), starting with u k and yielding the 
final result u k 

The cycle with 7 = 1 is called V cycle or V (17, 1/2), and the one with 7 = 2 
is called W cycle or W (17, 1/2) 


6.2 Full Multigrid (FMG) The JV-FMG is an algorithm for calculating an approx- 
imate solution 

u% = FMG{k,F k ,N) (6.4) 

to equation (6.1), defined recursively as the following two successive steps 

(a) Calculating a first approximation u k - If k = 1, put ufa = 0 Otherwise put 

u k = n\_ 1 FMG{k - l,l k k -'F k ,N), (6.5) 
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where 27* _! is an interpolation of solutions from grid h*_] to grid h*, and /* _1 is 
some transfer (“averaging”) from grid k to grid k — 1. usually a full weighting of 
the type (6 3). The interpolation 2 T£_j should usually be of order higher than that 
of the correction interpolation /*_j mentioned above [3, §7.1]. In our experiments 
bicubic interpolation was used 

(b) Improve the first approximation by N successive MG cycles 
tx? - MG(k,u k J _ i .F k ), (j = 1, . . N) 
as defined in Sec 6.1. 


6.3 Averaging The algorithm above is the conventional one, and for equations 
with constant coefficients it requires no modifications In case of quasi-elliptic 
equations with variable coefficients, and in particular in case of nonlinear equa- 
tions. it is not enough to prevent unstable initial errors, because such errors can 
also later be introduced due to interaction between modes It is then better to 
explicitly reduce the unstable modes by averaging, such as (3 1) or (3 2) It may 
also then be important to replace /*_ 1 u k ~ 1 in step (C) above by I k _ 1 S hl, - y u k ~ 1 
In fact, experiments with non-staggered Navier-Stokes equations (cf Sec 8 2) 
gave slowly diverging MG cycles unless this averaging was used 

6.4 Measuring convergence In various situations where algebraic convergence is 
not attempted, as in the present algorithm and double discretization [3, §10 2j 
and other algorithms, the question is raised how to measure convergence, how to 
know, m particular, that a solution to the truncation level (i e , with algebraic 
errors dominated by discretization errors) has been obtained 

The answer is that solution to the truncation level is not really the important 
information when differential convergence is our objective (as it should most often 
be), because (l) Solving to truncation level tells us nothing about the trunca- 
tion error itself We may for instance be doing good job in solving the algebraic 
system due to having chosen an easy-to-solve but badly-approximating discretiza- 
tion (n) A smaller differential error may often be obtained faster by switching to 
a finer grid before the equations on the present grid have been solved to truncation 
level 

The important information is the differential convergence itself as function 
of computational work This very information can directly be obtained from the 
TV-FMG algorithm Indeed, the sequence of approximations u k N , ( k = 1,2, ) is a 

sequence converging to the differential solution, hence the decrease in the sequence 
of differences 6* =jj 2* _1 ixjv — u^T 1 |j exactly exhibit the speed of differential 
convergence, where the norm || || used to measure <5* can be chosen to exactly 

represent the sense m which convergence is sought. One only has to check that 
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the smallness of 4 is not governed by lack of change from u k to It is enough 
for this purpose to check that the suitable residual norm r^- =]| F k — L k u k N || is 
considerably smaller than r£ One can usually also verify that the algebraic errors 
are below truncation level, e g , by confirming that /r k is considerably smaller 
than 


7 Algorithm for Fast Algebraic Convergence 


Although fast asymptotic algebraic convergence is not needed for fast differen- 
tial convergence, it can still be produced by a more involved multigrid algorithm 
This algorithm (also described in [6]) may be interesting in its own right, since it is 
the simplest example of a new kind of algorithms (first mentioned m (1, §3.2] , and 
more fully in [81) for solving problems w r ith highly-oscillatory solutions, including 
highly indefinite problems 


7 1 Multiple coarse grid corrections Let , 0 2 . ,0* be all the components for 

which L h vanishes, or, more generally, the centers of all neighborhoods in which 
L h (9) is small Usually 0 1 = 0 Then (by [3, §1 lj, for example) there exists a 
relaxation scheme with fast convergence for all Fourier components except those 
close to some 9 3 The error after few such relaxation sweeps must therefore have 
the form 

l 

^(i) = ]TV/*(i)exp {i9 J x/h ). (7.1) 

j=i 

where V k are smooth functions. Whereas classical multigrid seeks to approximate 
V h on a coarser grid and the algorithm of Sec 6 approximates Vj 1 , the new’ 
algorithm will separately approximate each of the V fc , by successively employing 
£ different coarse-grid corrections 

Generally, denoting by H the coarser-grid meshsize (H = 2 h), the equations 
for Vj H , the coaxse-grid approximation to Vj 1 , should have the form 

(7 2) 

where (0) & L k {9 3 + 0) for small 0, and //^(l?*) ~ 4* (=0 except for 6 }J = 1) 
The boundary conditions may couple and V^ 1 on any piece of boundary along 
which exp (t(9 3 —6 h ) x/h) is a smooth function There are various ways, variational 
ones and more direct ones to derive , lfj 3 and the boundary conditions There 
also exist various ways for solving (7.2). In highly indefinite problems the latter 
leads to creating more components on grid 4 h, etc., so that on increasingly coarser 
grids the representation tends to a Fourier representation 
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Here we give only the very simple example of solving for the skew Laplacian 
(2.5). (For a more general case, see [6] ) In this case £ = 2. 0 1 = (0,0) and 
6 2 = (7r, 7r ) , and one may simply take L? = to be any //-approximation to the 
Laplace operator In some situations, w’here the same mechanism creates both the 
fine grid and the coarse grid equations, these L ^ may again be skew-Laplacians. 
As transfer operators one can use 


I H - — 
lh 1 " 16 


1 2 1 

2 4 2 

.1 2 1J 


and — — 

h ’ 2 16 


1 -2 1 
-2 4 -2 

L 1 -2 1 


(7.3) 


Considering the case that the fine-grid boundary conditions are Dirichlet condi- 
tions identically satisfied by any fine-grid approximation, the coarse-grid boundary 
conditions for both Vf* and V 2 W are the homogeneous Dirichlet conditions For 
solving the coarse-grid equations (7 2) the MG cycle of Sec 6 1 can be used, even 
tn the case that L ^ are themselves quasi- elliptic, because, for the purpose of accel- 
erating the fine-grid algebraic convergence, equations (7 2) need to be solved each 
time only to their truncation level (l e , only to the level of the error V } H — V h ) 
In case of similar equations but w r ith non-constant coefficient, averaging as m Sec 
6 3 should better be used 


7.2 The modified algorithm Given an approximate solution to (6.1), the mod- 
ified multigrid cycle MMG for producing the improved solution u k 

u k <- MMG{k,u^F k ) (7 4) 

is defined nonr-recursively as follows 

If k = 1 solve (6 1) by any direct or iterative method, yielding the final u k 
Otherwise, perform u relaxation sweeps on (6.1), resulting m a new approximation 
u k '°, and then, for j = 1,2. . , £, calculate 

— MG{k- 1,0, /j^* [F k - L k u k ^~ 1 )) 

u k ' 3 = + exp(z0 J x/h) I k _ i v k ~ 1,3 

with u k = u k,e being the final result /£_ j again denotes linear interpolation MG 
is the cycle defined m Sec 6.1, with a choice of 7, U\. j/ 2 

With this MMG cycle replacing the MG cycle, the modified FMG algorithm 
is defined m the same way as FMG in Sec 6.2 

7.3 Modified smoothing analysis. The smoothing factor for the above MMG cycle, 
i.e., the ideal factor of convergence one can expect from such a cycle per relaxation 
sweep on the finest grid is defined by 

p = max l^(£)!- (7 5) 

ir/2<ie-e j | for one j, ie|<w 
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where |/i(0)| is the spectral radius of the amplification matrix (or the absolute value 

of the amplification factor, if q = 1) Note that for l > 2 d , the domain of 0 over 

which the maximum is taken may be empty. In such a situation convergence can 

m principle be obtained without any relaxation on the finest grid This does not 

mean that the algorithm is more efficient than a conventional multigrid, because 

it employs at least t times as many relaxation sweeps on each coarser grid 

# 

A more precise two-level analysis can of course be made here in the conven- 
tional way [3, §4 1] 

For the skew-Laplacian and the algorithm described above, the lexicographic 
Gauss-Seidel amplification factor (5 2) attains its maximum (7.5) at (±7t/2,0) and 
at (±7r/2,7r), yielding fi — 447 


8 Numerical Experiments 

8 1 The skew Laplacian problem Our main experimental studies were conducted 
with the skew Laplacian scheme (2.5) in the rectangle {0 < < 2,0 < 12 < 3} 

with Dirichlet boundary conditions These conditions and the right-hand side 
of the differential equation A U = F were chosen so that the solution U of the 
differential equations is known, to allow direct measurements of discretization 
errors The sequence of grids have meshsizes hk = 2 1 ~ k ( k = 1.2, .), each 

positioned so that the boundaries of fl coincide with grid lines On every level 
L k is the skew Laplacian, and the relaxation is lexicographic Gauss-Seidel The 
algorithms were those described m Secs 6 and 7. 

Table 1 shows the maximal differential error (maximal differences between 
computed and differential solutions) on various grids In addition, columns headed 
by d or d c show maximal error in first derivatives, approximated either at grid 
midpoint by short difference quotients (the d columns), or at gndpoints by d c (the 
d c columns). The upper part of the table gives these errors for the exact discrete 
solution, the lower part - for the solution obtained by a 1-FMG algorithm with 
V (2, 1) cycles For grid 5 an additional result (5a) is sometimes given It shows 
errors measured after the solution is averaged by ( T + T* 1/2 )/ 2 (cf. (3 1)) 
The table compares skew-Laplacian with usual (compact) Laplacian (using the 
same meshsize and the same relaxation), and a case of smooth solution with a 
highly-oscillatory case The latter is shown in order to emphasize how bad quasi- 
elliptic schemes can be. In practice such highly oscillatory components have very 
small amplitudes If their amplitudes are bigger than 0(h 2 ) (here h 2 — .001), 
then second-order approximations cannot be obtained by any discretization In 
the highly oscillating case it was of course necessary to use the full weighting (6 3) 
for m (6.5), this was started with k = 7. In the smooth case, however, 

injection of F was used, in order to obtain a clearer picture, clean of F-averaging 
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8 2 The Stokes and Navier- Stokes problems We have also conducted experiments 
with the Stokes operator (2.3), described in detail in [3, §18.6] (with slight improve- 
ments, to be described m the new edition) The unknown grid functions of this 
operator are U h , V h and P h - the discrete horizontal velocity, vertical velocity 
and pressure, respectively 

In the differential problem only velocities are normally given on the boundary. 
In the non-staggered discretization (2 3) some boundary conditions for P h should 
be introduced (which is a disadvantage typical to many quasi-elliptic operators) 
For clarity of exposition we here avoid this issue by showing results for periodic 
boundary conditions (adjusting undetermined additive constants before measuring 
errors). 

The exact treatment of boundary conditions is important only m measuring 
asymptotic convergence rates It does not much affect results of I-FMG Therefore 
we will show such results also for the Dirichlet boundary conditions In these 
experiments P h at each boundary point is taken equal to the nearest interior 
value of P h , and it changes whenever the latter does This does not correspond 
to Neumann boundary conditions, but to coupling the four subgrids into which 
the P h grid decouples. A partial relaxation sweep near Dirichlet boundaries is 
performed before each full relaxation sweep 

The relaxation employed is distributed Gauss-Seidel (DGS), a special case of 
a scheme for relaxing general PDE systems, explained in [3. §3.7] Briefly, it is 
equivalent to writing U h = - dltp%, V h = ~ ^ 2^3 and P h = — A^ 1 ^, 
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and relaxing by usual Gauss-Seidel the resulting equations in The changes 
m the latter imply changes mU h , V h and P h , which define the actual changes 
performed by the DGS relaxation The relaxation ordering is 4-colored, relaxing 
the four mentioned subgrids one at a time 

The domain for this problem is the square {0 < x 3 < 2n} The meshsizes are 
hk = 2~ fc 7r The right-hand side and the boundary' conditions are chosen so as to 
give the prescribed solution U — V = P = sin(cos(xi + 2 x 2 )), a periodic solution 
which includes many Fourier modes The discrete right-hand sides were calculated 
by F fc_1 = F k , using (6 3), starting at k = 8 

Some experiments were conducted with averaging (cf Sec 6 3). In the 
present case this means averaging of P h only, since U h and V h vanish m the 
unstable modes When used, this P-averagmg employed (3.1), with mj = m 2 = 2. 
performed on P h in any solution or correction just before interpolating it to a finer 
grid 

Also mentioned below are experiments with non-stagered incompressible 
Navier-Stokes ( INS) equations, with procedures similar to those for Stokes For 
details see [3. §19j the modification from staggered to non-staggered formulation 
and processing are the same as for the Stokes equations Results are given for the 
Dirichlet problem (U and V given on the boundary, P on the boundary treated as 
above), for the case U = V = P = l + .2sin(cos(x 1 +2x 2 )). We have experimented 
with small and large Reynolds numbers. Re In the latter case anisotropic artificial 
viscosity was used in relaxation, its magnitude being 1 4 times the viscosity intro- 
duced by upstream differencing Central differencing without artificial viscosity 
was used for the fine-to-coarse residual calculations, allowing 0(h 2 ) solutions to be 
obtained The large Re PDE problem is not elliptic (more precisely, it has small 
ellipticity measure), so its detailed discussion is beyond the scope here Indeed, 
the present example is not fully typical for large Re, because it has no boundary 
layers and no gridline-streamline alignments 

Table 2 summarizes four numerical experiments Three with the Stokes (Re = 
0) problem (exact solutions for the periodic (“Per ”) boundary conditions, 1-FMG 
solutions with W{ 2, 1) cycles for the same problem, and similar 1-FMG solutions 
for the Dirichlet (“Dir.”) problem), and one experiment for “infinite” Re. ie., 
with viscosity completely dominated by artificial viscosity The latter experiment 
uses 2-FMG algorithm with W (2, 0) cycles because double discretization (differ- 
ent artificial viscosities at different stages) is involved (cf. [3, §10.2]). For each 
experiment and each grid k , the three numbers shown in the first column are 
max(|| u k — U 1], |j v k — V ||), || p k - P || and j| p k — P ||, where ( u k ,v k ,p k ) is 
the solution obtained for that grid, p k = ^II^Lj (T^ 2 + T~ 1 ^ 2 )p k and || || is the 

discrete L\ norm per unit area. The three numbers in the next column (headed 
by “d”) are max J=lt2 max(|| d k u k - d 3 U ||, || d k v k - d 3 V ||), || d k p k - d 3 P ||, and 

j| d k p k - djP ||, where d 3 = d/dx 3 and d k = {T^ 2 — T~ l ^ 2 )/hk In the next col- 
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Re = 0 | 5 | 00084 0030 0108 .00084 .0030 0108 00079 .0055 
Per j | 00500 .0394 0392 .00500 0394 0392 00150 .0094 

I I I 


00542 .0395 .0395 .00542 0395 0395 00513 0212 

! Exact j 6 00024 .0007 0026 j 00024 0007 0026 00018 .0013 

1 I 

; Sol j I 00123 0098 0098 1 00123 0098 .0098 00037 0025 i 

I i I 


! 00135 0100 0108! 00135 0100 0108 j 00127 0054 
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umn (headed by d c ). similar numbers are given, with the long difference quotient 
d k — d' = (T, — T~ i )/(2hk) replacing d k The next 3 columns show similar sets 
of results for the case that P-averagmg is used The remaining 2 columns give for 
comparison results obtained on a staggered grid with the same meshsize, without 
P-averaging (Using p k for approximation, especially of derivatives, still may pay 
if 1^2 = 0) 


8 3 Accuracy and Stability Tables 1 and 2 clearly show that the exact quasi- 
elliptic solutions (A v and the non-staggered Stokes, the latter mainly in terms 
of P) are several times less accurate than the corresponding fully elliptic ones 
(A h and staggered Stokes, respectively), but they are still 0(h 2 ) Errors in the 
highly oscillating case exhibiting instability, could of course all be reduced to 0(1) 
(or 0(h7 1 ) in derivatives) by enough F-averagmg (see §3) Averaging the solu- 
tion (row 5a, or the p k results), or taking suitable long difference quotients, cure 
the worst behavior too but also somewhat further reduce the smooth-component 
accuracy, which nevertheless remains 0(h?) 


8 4 Poor asymptotic algebraic convergence Denote by A the asvmptotic conver- 
gence factor per multigrid cycle, l e . A = (r£/r m ) 1/, ^~ m ^ for sufficiently large £, 
m and i - m, where t-£ is any error (or residual) norm measured at any fixed stage 
of the f-th cycle As expected (see §4 1). the usual cycles MG(k ) yielded poor 
X for quasi-elliptic schemes 

In case of the Skew Laplacian and V (2 1) cycles, our experiments exhibited 
X = .845 and X = 96 for levels k = 4 and k = 5, respectively The convergence 
rate log 1/A is clearly 0(h 2 ), as the rate of a simple Gauss-Seidel solver for the 
compact Laplacian A h Indeed, on each subgrid (red or black) the relaxation does 
look like Gauss-Seidel for A hl . and the coarse grid corrections are no help m case 
the black residuals cancel the red ones in the transfer to grid k— 1 For comparison 
V(2, 1) cycles for the compact Laplacian A h with lexicographic Gauss-Seidel yield 
A % 12 on all grids 

Similarly, for the periodic Stokes problem and W(2. 1) cycles, A = 80 and 
A = 945 were obtained on levels 4 and 5, respectively, exhibiting again 0(h 2 ) 
rate. The rates were almost identically the same whether P averaging was used 
or not. For comparison, for staggered- grid Stokes discretizations the red-black 
DGS relaxation gives A = 30 and A = 20 for the IF(1,0) and the W'(2,0) cycles, 
respectively These same excellent rates are obtained both for the periodic and 
the Dirichlet boundary conditions (provided some local relaxation near boundaries 
is added in the latter case) The same results are obtained for the Navier-Stokes 
problem with small Re For large Re, divergence occur unless P-averaging is used 
(cf §6 3). 
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8 5 FMG results Despite the bad asymptotic convergence, Tables 1 and 2 clearly 
show that results obtained for the quasi-elhptic cases by short FMG algorithms 
are very good In smooth cases they yield differential errors practically as small 
as in the exact discrete solutions Moreover, in case of the unstable mode, the 
FMG results are visibly much better than the exact solution (precisely because the 
bad behavior is slow to enter) In case of non-linear equations (Table 2, Re — oc ) 
proper averaging (Sec 6 3) is Evidently necessary for good FMG results 


8 6 Asvmptotic comergence with new algorithm The MMG{b ) cycle of §7 2 
has been employed to solve the skew Laplacian problem with v — 3 relaxation 
sweeps per cycle and with V'(2, 1) used as the MG [A. I inner cycle For many 

cycles the convergence factor per cycle was steadily' between 07 and 08, or a 
convergence factor of 425 per fine-grid relaxation close to the value 447 expected 
by the smoothing mode analysis (§7 3) 
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